# Send message to log
log_code()

# figures ####
load("out/full2.Rda")
tmp_df <- full_df %>% 
  filter(parameter == "theta") %>% 
  select(mean, party, p_id, variable, country) %>% 
  # mutate(dim = rep(c("d1", "d2"), length.out = 256)) %>%
  mutate(dim = paste0("d", substr(variable, nchar(variable) - 1, nchar(variable) - 1))) %>%
  select(-variable) %>% 
  pivot_wider(names_from = dim, values_from = mean) %>% 
  merge(p_start)

cor(tmp_df$d1, tmp_df$lrecon)
cor(tmp_df$d2, tmp_df$galtan)

# figure 1 in main
left_join(tmp_df %>% 
            select(p_id, party, country, d1, d2) %>% 
            pivot_longer(cols = c(d1, d2)) %>% 
            rename(latent = value) %>% 
            mutate(name = ifelse(name == "d1", "Economic Dimension", "GAL-TAN Dimension")),
          tmp_df %>% 
            select(p_id, party, country, lrecon, galtan) %>% 
            pivot_longer(cols = c(lrecon, galtan)) %>% 
            rename(experts = value) %>% 
            mutate(name = ifelse(name == "lrecon", "Economic Dimension", "GAL-TAN Dimension"))) %>% 
  ggplot(aes(experts, latent)) + geom_point() + geom_smooth(method = "lm") + theme_bw() + geom_abline(aes(intercept = 0, slope = 1)) + facet_wrap(~name) +
  xlab("Mean Expert Placement") + ylab("Latent Estimate")
# ggsave(file = "figures/cor_ideal.pdf", width = 3*1.61, height = 3)
ggsave(file = "figures/fig1.pdf", width = 3*1.61, height = 3)

# figure 2 in main and figures 4, 5 in si
tmp_df <- full_df %>% 
  filter(parameter %in% c("beta_1", "beta_2")) %>% 
  mutate(item.dimension = ifelse(item %in% c("eu_asylum", "eu_budgets", "eu_cohesion", "eu_foreign", "eu_intmark", "eu_position"), "EU Issues",
                                 ifelse(item %in% c("civlib_laworder", "environment", "ethnic_minorities", "immigrate_policy", "multiculturalism", "nationalism", "regions", "religious_principles", "sociallifestyle", "urban_rural"), 
                                        "Social/Cultural Issues", "Economic Issues")))

filter(tmp_df, parameter %in% c("beta_1", "beta_2")) %>% 
  mutate(parameter = ifelse(parameter == "beta_1", "Economic Dimension", "GAL-TAN Dimension"),
         item.dimension = factor(item.dimension, levels = c("Economic Issues", "Social/Cultural Issues", "EU Issues"))) %>%
  ggplot(aes((item), mean)) + geom_hline(aes(yintercept = 0), linetype = 2) + geom_boxplot() + facet_grid(item.dimension~parameter, scales = "free", space = "free") + 
  coord_flip() + theme_bw() + theme(legend.position = "") + ylab("Discrimination Estimate") + xlab("Issue") + labs(color = "")
# ggsave("figures/disc_box.pdf", width = 5*1.61, height = 5)
ggsave("figures/fig2.pdf", width = 5*1.61, height = 5)

filter(tmp_df, parameter %in% c("beta_1", "beta_2")) %>% 
  mutate(parameter = ifelse(parameter == "beta_1", "Economic Dimension", "GAL-TAN Dimension"),
         item.dimension = factor(item.dimension, levels = c("Economic Issues", "Social/Cultural Issues", "EU Issues"))) %>% 
  filter(parameter == "Economic Dimension") %>% 
  ggplot(aes(info, mean, ymin = q5, ymax = q95)) + geom_linerange() + geom_point() + theme_bw() + facet_wrap(~item, ncol = 4) + xlab("Country") + ylab("Discrimination Estimate") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + geom_hline(aes(yintercept = 0))
# ggsave("figures/disc_econ.pdf", width = 4*1.61, height = 7)
ggsave("figures/fig4.pdf", width = 4*1.61, height = 7)

filter(tmp_df, parameter %in% c("beta_1", "beta_2")) %>% 
  mutate(parameter = ifelse(parameter == "beta_1", "Economic Dimension", "GAL-TAN Dimension"),
         item.dimension = factor(item.dimension, levels = c("Economic Issues", "Social/Cultural Issues", "EU Issues"))) %>% 
  filter(parameter != "Economic Dimension") %>% 
  ggplot(aes(info, mean, ymin = q5, ymax = q95)) + geom_linerange() + geom_point() + theme_bw() + facet_wrap(~item, ncol = 4) + xlab("Country") + ylab("Discrimination Estimate") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + geom_hline(aes(yintercept = 0))
# ggsave("figures/disc_galtan.pdf", width = 4*1.61, height = 7)
ggsave("figures/fig5.pdf", width = 4*1.61, height = 7)

# figure 3 in main
filter(full_df, parameter %in% c("idio")) %>% 
  mutate(item.dimension = ifelse(item %in% c("eu_asylum", "eu_budgets", "eu_cohesion", "eu_foreign", "eu_intmark", "eu_position"), "EU Issues",
                                 ifelse(item %in% c("civlib_laworder", "environment", "ethnic_minorities", "immigrate_policy", "multiculturalism", "nationalism", "regions", "religious_principles", "sociallifestyle", "urban_rural"), 
                                        "Social/Cultural Issues", "Economic Issues"))) %>% 
  mutate(item = fct_reorder(item, mean, .fun='sd')) %>% 
  ggplot(aes(item, mean, fill = item.dimension)) + geom_hline(aes(yintercept = 0), linetype = 2) + geom_boxplot() + 
  coord_flip() + theme_bw() + theme(legend.position = "bottom") + ylab("Idiosyncratic Shock") + xlab("Issue") + 
  scale_fill_grey(start = .4, end = .9) + labs(fill = "")
# ggsave("figures/idio.pdf", width = 4*1.61, height = 4)
ggsave("figures/fig3.pdf", width = 4*1.61, height = 4)

# figure 6 in si
filter(full_df,parameter %in% c("idio")) %>% 
  mutate(item.dimension = ifelse(item %in% c("eu_asylum", "eu_budgets", "eu_cohesion", "eu_foreign", "eu_intmark", "eu_position"), "EU Issues",
                                 ifelse(item %in% c("civlib_laworder", "environment", "ethnic_minorities", "immigrate_policy", "multiculturalism", "nationalism", "regions", "religious_principles", "sociallifestyle", "urban_rural"), 
                                        "Social/Cultural Issues", "Economic Issues"))) %>% 
  mutate(item = fct_reorder(item, mean, .fun='sd')) %>% 
  merge(df_pid %>% select(party = party_id, family)) %>% 
  merge(data.frame(family = 1:11, family.name = c("Radical Right", "Conservatives", "Liberal", "Christian-Democratic", "Socialist", "Radical Left", "Green", "Regionalist", "No Family", "Confessional", "Agrarian/Center"))) %>% 
  filter(family %in% c(1:8, 10:11)) %>% 
  ggplot(aes(item, mean, fill = item.dimension)) + geom_hline(aes(yintercept = 0), linetype = 2) + geom_boxplot() + facet_wrap(~family.name, ncol = 3) +
  coord_flip() + theme_bw() + theme(legend.position = "bottom") + ylab("Mean Idiosyncratic Shock") + xlab("Issue") + 
  scale_fill_grey(start = .4, end = .9) + labs(fill = "")
# ggsave("figures/idio_party.pdf", width = 9, height = 12)
ggsave("figures/fig6.pdf", width = 9, height = 12)

# figures 7:10 in si
full_df %>% 
  filter(parameter == "res_alpha") %>% 
  select(mean, q5, q95, country, item) %>% 
  mutate(item = fct_reorder(item, mean, .fun='sd')) %>% 
  mutate(item.dimension = ifelse(item %in% c("eu_asylum", "eu_budgets", "eu_cohesion", "eu_foreign", "eu_intmark", "eu_position"), "EU Issues",
                                 ifelse(item %in% c("civlib_laworder", "environment", "ethnic_minorities", "immigrate_policy", "multiculturalism", "nationalism", "regions", "religious_principles", "sociallifestyle", "urban_rural"), 
                                        "Social/Cultural Issues", "Economic Issues"))) %>% 
  ggplot(aes(item, mean, ymin = q5, ymax = q95, fill = item.dimension)) + geom_hline(yintercept = 0) + geom_boxplot() + theme_bw() + coord_flip() +
  xlab("Issue") + ylab("Estimate") + theme(legend.position = "bottom") +
  scale_fill_grey(start = .4, end = .9) + labs(fill = "")
# ggsave(file = "figures/shift_issue_country.pdf", width = 4*1.61, height = 4)
ggsave(file = "figures/fig7.pdf", width = 4*1.61, height = 4)

full_df %>% 
  filter(parameter == "scale_beta") %>% 
  select(mean, q5, q95, info) %>% 
  mutate(info = fct_reorder(info, mean, .fun='mean')) %>% 
  ggplot(aes(info, mean, ymin = q5, ymax = q95)) + geom_pointrange() + theme_bw() + 
  geom_hline(yintercept = 1) + xlab("Country") + ylab("Estimate")
# ggsave(file = "figures/scale_country.pdf", width = 2.25*1.61, height = 2.25)
ggsave(file = "figures/fig8.pdf", width = 2.25*1.61, height = 2.25)

full_df %>% 
  filter(parameter == "shift_alpha") %>% 
  select(mean, q5, q95, info) %>% 
  mutate(info = fct_reorder(info, mean, .fun='mean')) %>% 
  ggplot(aes(info, mean, ymin = q5, ymax = q95)) + geom_pointrange() + theme_bw() + 
  geom_hline(yintercept = 0) + xlab("Country") + ylab("Estimate")
# ggsave(file = "figures/shift_country.pdf", width = 2.25*1.61, height = 2.25)
ggsave(file = "figures/fig9.pdf", width = 2.25*1.61, height = 2.25)

full_df %>% 
  filter(parameter == "exp_scale") %>% 
  separate(info, c("country", "expert")) %>% 
  select(mean, country, expert) %>% 
  mutate(country = fct_reorder(country, mean, .fun='median')) %>% 
  ggplot(aes(country, mean)) + geom_hline(yintercept = 1) + geom_boxplot() + theme_bw() + xlab("Country") + ylab("Estimate")
# ggsave(file = "figures/expert_error.pdf", width = 2.25*1.61, height = 2.25)
ggsave(file = "figures/fig10.pdf", width = 2.25*1.61, height = 2.25)

rm(list = ls())

# model comparison ####
load("out/comp.Rda")

change_df <- data.frame(old = paste0("m", 0:5),
                        new = paste0("Model ", c(1, 4, 5, 2, 3, 6)))

loo  %>% 
  select(old = model, Estimate, SE) %>% 
  merge(change_df) %>% 
  arrange(new) %>% 
  ggplot(aes(as.character(1:6), Estimate, ymin = Estimate - 1.96 * SE, ymax = Estimate + 1.96 * SE)) + geom_pointrange() + theme_bw() + xlab("Model") + ylab("ELPD-LOO")
# ggsave(file = "figures/loo.pdf", width = 3, height = 2.5)
ggsave(file = "figures/fig11b.pdf", width = 3, height = 2.5)

p_acc  %>% 
  select(old = model, mean, q5, q95) %>% 
  merge(change_df) %>% 
  arrange(new) %>% 
  ggplot(aes(as.character(1:6), mean, ymin = q5, ymax = q95)) + geom_pointrange() + theme_bw() + xlab("Model") + ylab("Predictive Accuracy")
# ggsave(file = "figures/p_acc.pdf", width = 3, height = 2.5)
ggsave(file = "figures/fig11a.pdf", width = 3, height = 2.5)

loo <- loo %>% 
  mutate(ELPD_LOO = paste0(round(Estimate, 3), " (", round(SE, 3), ")")) %>% 
  select(old = model, ELPD_LOO) %>% 
  merge(change_df) %>% 
  select(-old)

p_acc <- p_acc %>% 
  mutate(Acc. = paste0(round(mean, 3), " (", round(q5, 3), "--", round(q95, 3), ")")) %>% 
  select(old = model, Acc.) %>% 
  merge(change_df) %>% 
  select(-old)

m_comp <- merge(loo, p_acc) %>% 
  rename(Model = new)
names(m_comp) <- c("Model", "ELPD-LOO", "Accuracy")

fileConn <- file("tables/mod_comp.tex")
writeLines(print(xtable(m_comp, 
                        align = c("cc|cc"),
                        caption = "Model comparison results",
                        label = "tab:mod_comp_res"), booktabs = TRUE, include.rownames = FALSE), fileConn)
close(fileConn)

